##-----------------------------------------##
## In-group similarity and social distance ##   
##-----------------------------------------##
rm(list = ls())


## Load data
yearly_between_group_distance <- readRDS("data/yearly_between_group_distance_replication.rds")
polity.diffusion <- readRDS("data/polity_diffusion_replication.rds")
war.data <- readRDS("data/war_data_replication.rds")

## Set up data to find polarization
predict.democ <- polity.diffusion
predict.democ$gov <- "Mixed regime"
predict.democ$gov[predict.democ$polity >= 6] <- "Democracy"
predict.democ$gov[predict.democ$polity <= -6] <- "Autocracy"
unique.groups <- unique(predict.democ$membership)

war.yearly <- plyr::ddply(
    war.data,
    ~year,
    summarize,
    mids = log(sum(fatal_mid) + 1)
  ) %>% 
  arrange(
    year
  ) 

groups.yearly <- distinct(
    polity.diffusion,
    membership, year
  ) %>% 
  mutate(
    n = 1
  ) %>% 
  plyr::ddply(
    ~year,
    summarize,
    groups = sum(n)
  )

states.yearly <- distinct(
    polity.diffusion,
    names, year
  ) %>% mutate(
    n = 1
  ) %>%
  plyr::ddply(
    ~year,
    summarize,
    states = sum(n)
  )

system_data <- filter(
    yearly_between_group_distance,
    year <= 2014
  ) %>% 
  inner_join(
    war.yearly,
    by = "year"
  ) %>% 
  inner_join(
    groups.yearly,
    by = "year"
  ) %>% 
  inner_join(
    states.yearly,
    by = "year"
  ) %>%
  mutate(
    log.groups = scale(log(groups + 1))
  )


polarization <- ts(data  = system_data$assortativity_regime_normalized,
                   start = 1816,
                   end   = 2014)
ingroup_ties <- ts(data  = system_data$assortativity_group_normalized, 
                   start = 1816,
                   end   = 2014)

groups <- ts(data  = system_data$groups, 
             start = 1816,
             end   = 2014)

states <- ts(data  = system_data$states, 
             start = 1816,
             end   = 2014)


prop_democ <- ts(data  = system_data$prop_democ, 
                 start = 1816,
                 end   = 2014)

modularity <- ts(data  = system_data$modularity, 
                 start = 1816,
                 end   = 2014)

mids <- ts(data  = log(war.yearly$mids + 1), 
           start = 1816,
           end   = 2014)


varDat <- data.frame(modularity, polarization, mids, ingroup_ties, groups, prop_democ, states)
varDat$outliers <- ifelse((1:nrow(varDat) + 1815) %in% c(1914:1918, 1939:1945, 1990), 1, 0)
varDat <- varDat[complete.cases(varDat),]
colnames(varDat) <- c("ingroup", "assortatative_regime", "mids", "assortatative_group", "groups", "prop_democ", "states", "systemic_wars")
exog <- data.frame(varDat[1:(nrow(varDat)),c(7)])
names(exog) <- "n_states"
varDatinput <- varDat[1:(nrow(varDat)),c(1:2)]
# varDatinput$ingroup <- as.numeric(scale(varDatinput$ingroup))
# varDatinput$assortatative_regime <- as.numeric(scale(varDatinput$assortatative_regime))

lagselection  <- VARselect(varDatinput, lag.max = 25, type = "trend", exogen = exog)
lagselection$selection

var1  <- VAR(varDatinput, p = 5,  type = "trend", exogen = exog)
set.seed(1)
granger_ingroup      <- causality(var1, cause = "ingroup", boot = T, boot.runs = 1000)
granger_ingroup$Granger
granger_assortatative_regime      <- causality(var1, cause = "assortatative_regime", boot = T, boot.runs = 1000)
granger_assortatative_regime$Granger
summary(var1)
save(granger_ingroup, file = "model_output/granger_ingroup.RData")
save(granger_assortatative_regime, file = "model_output/granger_assortatative_regime.RData")
save(var1, file = "model_output/var_assortativity.RData")


granger_in_group <- c()
granger_assort <- c()
for (i in 1:10)
{
  var1  <- VAR(varDatinput, p = i,  type = "trend", exogen = exog)
  set.seed(1)
  granger_ingroup      <- causality(var1, cause = "ingroup", boot = T, boot.runs = 1000)
  granger_in_group[i] <- granger_ingroup$Granger$p.value
  granger_assortatative_regime      <- causality(var1, cause = "assortatative_regime", boot = T, boot.runs = 1000)
  granger_assort[i] <- granger_assortatative_regime$Granger$p.value
  save(granger_ingroup, file = paste0("model_output/granger_ingroup_lag_", i, ".RData"))
  save(granger_assortatative_regime, file = paste0("model_output/granger_assortatative_regime_lag_", i, ".RData"))
  save(var1, file = paste0("model_output/var_assortativity_lag_", i, ".RData"))  
}



yearly_between_group_distance <- arrange(
  yearly_between_group_distance, 
  year
)
yearly_between_group_distance$mids <- mids
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
plot_assortativity <- mutate(
    yearly_between_group_distance,
    modularity = range01(modularity), 
    assortativity_regime = range01(assortativity_regime)
  ) %>% 
  reshape2::melt( 
    id = c("year")
  ) %>% 
  filter(
    variable %in% c("modularity", "assortativity_regime")
  ) %>% 
  mutate(
    variable = as.character(variable),
    variable = ifelse(grepl("regime", variable), "Same\nregime ties", variable),
    variable = ifelse(grepl("mids", variable), "Fatal MIDs\nover states", variable),
    variable = ifelse(grepl("modular", variable), "Modularity", variable)
  ) %>% 
  mutate(
    title = "Modularity and same regime assortativity"
  )
saveRDS(plot_assortativity, file = "data/plot_assortativity.rds")

